home *** CD-ROM | disk | FTP | other *** search
/ NeXT Education Software Sampler 1992 Fall / NeXT Education Software Sampler 1992 Fall.iso / Mathematics / Notebooks / SigProc2.0 / Packages / SignalProcessing / Digital / DTFT.m < prev    next >
Encoding:
Text File  |  1992-08-18  |  36.4 KB  |  1,122 lines
  1. (*  :Title:    Discrete-Time Fourier Transforms *)
  2.  
  3. (*  :Authors:    Wally McClure, Brian Evans, James McClellan  *)
  4.  
  5. (*
  6.     :Summary:    This file will compute the forward and inverse
  7.         discrete-time Fourier Transform of a signal.
  8.  *)
  9.  
  10. (*  :Context:    SignalProcessing`Digital`DTFT`     *)
  11.  
  12. (*  :PackageVersion:  2.7    *)
  13.  
  14. (*
  15.     :Copyright:    Copyright 1989-1991 by Brian L. Evans
  16.         Georgia Tech Research Corporation
  17.  
  18.     Permission to use, copy, modify, and distribute this software
  19.     and its documentation for any purpose and without fee is
  20.     hereby granted, provided that the above copyright notice
  21.     appear in all copies and that both that copyright notice and
  22.     this permission notice appear in supporting documentation,
  23.     and that the name of the Georgia Tech Research Corporation,
  24.     Georgia Tech, or Georgia Institute of Technology not be used
  25.     in advertising or publicity pertaining to distribution of the
  26.     software without specific, written prior permission.  Georgia
  27.     Tech makes no representations about the suitability of this
  28.     software for any purpose.  It is provided "as is" without
  29.     express or implied warranty.
  30.  *)
  31.  
  32. (*
  33.     :History:    start            October 2, 1989
  34.         made into package    April 18, 1990
  35.         made into release    May 25, 1990
  36.         debugged fully        September 6, 1990
  37.  *)
  38.  
  39. (*  :Keywords:    discrete-time Fourier transform  *)
  40.  
  41. (*
  42.     :Source:    {Discrete-Time Signal Processing} by A. V. Oppenheim and
  43.           Ronald W. Schafer.  Prentice Hall.  1989
  44.  *)
  45.  
  46. (*
  47.     :Warning:    If the input does not match any given transform, a nasty
  48.         expression is returned. (But isn't this true for any
  49.         invalid transform in the signal processing packages?)
  50.  *)
  51.  
  52. (*  :Mathematica Version:  1.2 or 2.0  *)
  53.  
  54. (*  :Limitation:  *)
  55.  
  56. (*
  57.     :Discussion:  This  file  implements the forward and inverse DTFT.
  58.           The DTFTransform rule base (DTFTRules) consists of 5
  59.           property rules, 2  structure  rules  (upsampling and
  60.           downsampling),  13 lookup rules (for those functions
  61.           which  have  a  DTFT  but  not a z-transform), and 5
  62.           strategies --   distribute polynomials, expand terms
  63.           in cosine, sine,  and  exponential functions,  and a
  64.           call to  the  forward  z-transform rule base (if all
  65.           else fails).  The only functions supporting the for-
  66.           ward DTFT which depend  on the implementation of the
  67.           z-transform are ZtoDTFT and DTFTtoZ.
  68.  
  69.           The InvDTFTransform is similarly arranged, having 11
  70.           lookup  rules,  6  property  rules, and 6 strategies
  71.           (including  a  call  to the inverse z-transform rule
  72.           base).  Only one rule in  InvDTFTRules  is dependent
  73.           on  the  implementation  of  the inverse z-transform
  74.           rule base InvZTransform.
  75.  
  76.     MyDTFT[f, n, w] determines the discrete-time Fourier transform
  77.     of f with respect to n, where w is the frequency variable.
  78.  
  79.     MyInvDTFT[f, w, n] determines the inverse discrete-time Fourier
  80.     transform of f with respect to w, where n is a "time" variable.
  81.  
  82.     These transform pairs are valid with respect to the DTFT formulas:
  83.  
  84.         n
  85.     (-1)                 w
  86.     -----  <----->  - ------ CPulse[2 Pi, w + Pi]
  87.       n               2 I Pi
  88.  
  89.               n
  90.               (-1)
  91.     - 2 I Pi  -----  <----->  w CPulse[2 Pi, w + Pi]
  92.                 n
  93.  *)
  94.  
  95. (*
  96.     :Functions:    DTFTransform
  97.         InvDTFTransform
  98.  *)
  99.  
  100.  
  101.  
  102. (*  B E G I N     P A C K A G E  *)
  103.  
  104. BeginPackage[ "SignalProcessing`Digital`DTFT`",
  105.           "SignalProcessing`Digital`DSupport`",
  106.           "SignalProcessing`Digital`InvZTransform`",
  107.           "SignalProcessing`Digital`ZTransform`",
  108.           "SignalProcessing`Digital`ZSupport`",
  109.           "SignalProcessing`Support`TransSupport`",
  110.           "SignalProcessing`Support`ROC`",
  111.           "SignalProcessing`Support`SigProc`",
  112.           "SignalProcessing`Support`SupCode`" ]
  113.  
  114.  
  115. If [ TrueQ[ $VersionNumber >= 2.0 ],
  116.      Off[ General::spell ];
  117.      Off[ General::spell1 ] ];
  118.  
  119.  
  120. (*  U S A G E     I N F O R M A T I O N  *)
  121.  
  122. DTFTData::usage =
  123.     "DTFTData is a data-tag."
  124.  
  125. DTFTransform::usage =
  126.     "DTFTransform[function, in variable(s), out variable(s), options] \
  127.     returns the discrete-time Fourier transform of the input function. \
  128.     Note that this transform may rely on the z-transform and the \
  129.     substitution z = Exp[I w], which is not valid if the region of \
  130.     convergence does not include the unit circle.  Also note that \
  131.     only the first two arguments are required."
  132.  
  133. InvDTFTransform::usage =
  134.     "InvDTFTransform[function, in variable(s), out variable(s), options] \
  135.     return the inverse discrete-time Fourier transform (inverse DTFT) \
  136.     of the input function. \
  137.     The first two arguments are required. \
  138.     By default, the inverse DTFT rule base will apply the definition \
  139.     of the inverse DTFT when all other attempts to find the inverse \
  140.     have failed."
  141.  
  142. (*  E N D     U S A G E     I N F O R M A T I O N  *)
  143.  
  144.  
  145. Begin[ "`Private`" ]
  146.  
  147.  
  148. (*  B E G I N     P A C K A G E  *)
  149.  
  150. (*  Global variables  *)
  151. nList = {}
  152. wList = {}
  153. zList = {}
  154. multidFlag = False
  155.  
  156. (*  displayFixUp  *)
  157. displayFixUp[ fixUp_, op_ ] :=
  158.     Append[ fixUp[op], TransformLookup -> Replace[TransformLookup, op] ]
  159.  
  160. (*  validvarQ-- Valid variable checker  *)
  161. validvarQ[ x_Symbol ] := True
  162. validvarQ[ x_List ] := Apply[And, Map[VariableQ, x]]
  163. validvarQ[ x_ ] := False
  164.  
  165. (*  Generate a unique symbol with a Global context.  *)
  166. unique[x_] :=
  167.     Block [    {context, newsymbol},
  168.         context = $Context;
  169.         $Context = "Global`";
  170.         newsymbol = Unique[x];
  171.         $Context = context;
  172.         newsymbol ]
  173.  
  174.  
  175. (*  B E G I N     F O R W A R D     D T F T  *)
  176.  
  177. (*  Conversion from operator notation to function notation  *)
  178. Unprotect[DTFT]
  179. DTFT/: TheFunction [ DTFT[n_, w_][f_] ] := DTFTransform[f, n, w]
  180. Protect[DTFT]
  181.  
  182. (*  Options for the DTFT  *)
  183. DTFTransform/: Options[ DTFTransform ] :=
  184.     displayFixUp[ fixUp, { Dialogue -> False } ~Join~ Options[ZTransform] ]
  185.  
  186. (*  Extension of TheFunction operator.  *)
  187. DTFTData/: TheFunction [ DTFTData[ trans_, var_ ] ] := trans
  188.  
  189. (*  Magnitude/Phase plot of a DTFT object *)
  190. mppoptions = {    Domain -> Continuous, DomainScale -> Linear,
  191.         MagRangeScale -> Linear, PhaseRangeScale -> Degree,
  192.         PlotRange -> All }
  193.  
  194. DTFTData/: MagPhasePlot[DTFTData[trans_, FVariables[w_Symbol]]] :=
  195.     MagPhasePlot[ trans, {w, -Pi, Pi}, mppoptions ]
  196.  
  197. DTFTData/: MagPhasePlot[DTFTData[trans_, FVariables[{w1_Symbol, w2_Symbol}]]] :=
  198.     MagPhasePlot[ trans, {w1, -Pi, Pi}, {w2, -Pi, Pi}, mppoptions ]
  199.  
  200. (*  Input a function and test to see if needs to be sent to    *)
  201. (*  the multidimensional or one-dimensional Fourier transform.    *)
  202. DTFTransform[ f_ ] :=
  203.     Message[ Transform::novariables, "n (time)", GetVariables[f] ]
  204.  
  205. DTFTransform[ f_, n_ ] :=
  206.     DTFTransform [ f, n, DummyVariables[Length[n], Global`w] ]
  207.  
  208. DTFTransform[ f_, n_, w_, options___ ] :=
  209.     Message[ Transform::badvar, "discrete-time", DTFTransform, n ] /;
  210.     ! validvarQ[n]
  211.  
  212. DTFTransform[ f_, n_, w_, options___ ] :=
  213.     Message[ Transform::badvar, "continuous-frequency", DTFTransform, w ] /;
  214.     ! validvarQ[w]
  215.  
  216. DTFTransform[ f1_, n_, w_, options___ ] :=
  217.     Block [    {op, trans, vars, zvars},
  218.         vars = If [ Length[n] != Length[w],
  219.                 DummyVariables[Length[n], Global`w],
  220.                 w ];
  221.         zvars = DummyVariables[Length[n], z];
  222.  
  223.         nList = ToList[n];
  224.         wList = ToList[vars];
  225.         zList = ToList[zvars];
  226.  
  227.         op = ToList[options] ~Join~ Options[DTFTransform];
  228.         multidFlag = ( Length[n] > 1 );
  229.         trans = If [ multidFlag, 
  230.                  multiDDTFT[f1, n, vars, zvars, op],
  231.                  MyDTFT[f1, n, vars, zvars, op] ];
  232.  
  233.         If [ InformUserQ[op],
  234.              Scan[explain, trans, Infinity] ];
  235.  
  236.         DTFTData[ trans, FVariables[vars] ] ]
  237.  
  238. (* Multidimensional discrete-time Fourier transform *)
  239. multiDDTFT[ f1_, nl_, wl_, zl_, oplist_ ]:= 
  240.     Block [    {fun = f1, i, len, newop},
  241.         len = Length[nl];
  242.         newop = MultiDLookup[oplist, nl, wl];
  243.         For [ i = 1, i <= len, i++,
  244.               fun = MyDTFT[fun, nl[[i]], wl[[i]], zl[[i]], newop] ];
  245.         fun ]
  246.  
  247. (*  One-dimensional discrete-time Fourier transform  *)
  248. MyDTFT[ x_, n_, w_, z_, oplist_ ] :=
  249.     Block [ {newexpr = myDTFT[x, n, w, z, fixUp[oplist], initDTFTstate[]],
  250.          newdtftrules, newtrans, oldexpr = Null, ret, trace},
  251.  
  252.         (* generate the DTFT rules *)
  253.         newdtftrules = TransformFixUp[ n, w, oplist, myDTFT, False,
  254.                            DTFTransform, Null, Null ];
  255.         DTFTRules = Join[ newdtftrules, OriginalDTFTRules ];
  256.  
  257.         (* determine the level of dialogue *)
  258.         trace = SameQ[ Replace[Dialogue, oplist], All ];
  259.  
  260.         (* take the 1-D DTFT *)
  261.         Off[ Power::infy ];
  262.         While [ ! SameQ[newexpr, oldexpr], 
  263.             If [ trace, Print[newexpr]; Print[ "which becomes" ] ];
  264.             oldexpr = newexpr;
  265.             newexpr = oldexpr /. ( myDTFT[a__] :>
  266.                     Replace[ myDTFT[a], DTFTRules ] ) ];
  267.  
  268.         (* fix up the result *)
  269.         newexpr = posttransform[oldexpr];
  270.         While [ ! SameQ[newexpr, oldexpr], 
  271.             If [ trace, Print[newexpr]; Print[ "which becomes" ] ];
  272.             oldexpr = newexpr;
  273.             newexpr = posttransform[oldexpr] ];
  274.  
  275.         (* convert any remaining z-transform artifacts  *)
  276.         If [ ! FreeQ[{newexpr}, z],
  277.              newexpr = ( newexpr /. z -> Exp[I w] ) ];
  278.         On[ Power::infy ];
  279.         If [ trace, Print[newexpr] ];
  280.         newtrans = If [ SameQ[Head[newexpr], ZTransData],
  281.                 TheFunction[newexpr],
  282.                 newexpr ];
  283.  
  284.         (* properly handle Periodic operators *)
  285.         newtrans = newtrans /.
  286.             Periodic[k_, e_][f_] :> Periodic[k, Expand[e]][f];
  287.         newtrans = newtrans //.
  288.             { a_ Periodic[k_, w][f_] :> Periodic[k, w][a f],
  289.               Periodic[k_, w + w0_][f_] :> Periodic[k, w][f] /;
  290.                 FreeQ[w0, w],        (* ignore shifts *)
  291.               Periodic[k_, a_ w + w0_.][f_] :>
  292.                 Periodic[k / a, w][f] /;
  293.                 FreeQ[w0, w] };        (* scale axis *)
  294.         If [ trace, Print[ "which becomes" ]; Print[newexpr] ];
  295.  
  296.         (* simplify expression if requested *)
  297.         If [ Replace[Simplify, op],
  298.              newtrans = SPSimplify[newtrans,
  299.                        Dialogue -> trace,
  300.                        Variables -> wList] ];
  301.  
  302.         newtrans ]
  303.  
  304. (*  Supporting Routines  *)
  305.  
  306. diagonalResamplingMatrixQ[mat_] := False  /; ! MatrixQ[mat]
  307. diagonalResamplingMatrixQ[mat_] :=
  308.     Block [ {cond = False, diag, dims, i},
  309.         dims = Dimensions[mat];
  310.         If [ dims[[1]] == dims[[2]],
  311.              diag = Table[ mat[[i,i]], { i, 1, dims[[1]] } ];
  312.              cond = Apply[And, Map[isinteger, diag] ] &&
  313.                 SameQ[mat, DiagonalMatrix[diag]] ];
  314.         cond ]
  315.  
  316. explain[ myDTFT[ x_, n_, rest__ ] ] :=
  317.     Message[ Transform::incomplete,
  318.          "forward discrete-time Fourier transform", x, n ]
  319.  
  320. (*  fixUp --  removes redundancy and TransformLookup in options list  *)
  321. fixUp[ op_ ] := { Definition -> Replace[Definition, op],
  322.           Dialogue -> Replace[Dialogue, op] }
  323.  
  324. isinteger[x_] := Implies[ NumberQ[x], IntegerQ[x] ]
  325.  
  326. isintegerMatrix[ x_ ] :=
  327.     MatrixQ[x] && Apply[ SameQ, Dimensions[x] ] &&
  328.     Apply[ And, Map[isinteger, Flatten[x]] ]
  329.  
  330. mDresampleCheck[Downsample, n_, nvars_, nList_, m_] :=
  331.     mDresampleCheck[Upsample, n, nvars, nList, m]
  332.  
  333. mDresampleCheck[Upsample, n_, nvars_, nList_, l_] :=
  334.     Block [    {cond},
  335.         cond = SubsetQ[{n}, nvars, nList] && isintegerMatrix[l];
  336.         If [ cond, Needs[ "SignalProcessing`Digital`MDDTFT`" ] ];
  337.         cond ]
  338.  
  339. posttransform[x_] := ( x /. postrules )
  340.  
  341. postrules = {
  342.     (*    Apply the definition    *)
  343.     myDTFT[ x_, n_, w_, z_, options_, flag_ ] :>
  344.         Block [ {def, summhead},
  345.                 summhead = If [ TrueQ[$VersionNumber >= 2.0],
  346.                     Needs[ "Algebra`SymbolicSum`" ];
  347.                     Algebra`SymbolicSum,
  348.                     Needs[ "Algebra`GosperSum`" ];
  349.                     Algebra`GosperSum ];
  350.             def = summhead[x Exp[I w n], {n, -Infinity, Infinity}];
  351.             TransformDialogue[ Definition, options,
  352.                        summhead, x, def ] ] /;
  353.         Replace[Definition, options],
  354.  
  355.     derivative[x_, w_, 1] :>
  356.         I D[ posttransform[x], w ],
  357.     derivative[x_, w_, k_] :>
  358.         I^k D[ posttransform[x], {w, k} ],
  359.     downsample[x_, w_Symbol, m_Integer] :>
  360.         Block [ {k, newx},
  361.             newx = posttransform[x] /. w -> ((w + 2 Pi k) / m);
  362.             Sum[newx, {k, 0, Abs[m]-1}] / Abs[m] ],
  363.     downsample[x_, w_Symbol, m_] :>
  364.         Block [ {k, newx},
  365.             k = Unique["k"];
  366.             newx = posttransform[x] /. w -> ((w + 2 Pi k) / m);
  367.             Summation[k, 0, Abs[m]-1, 1][newx] / Abs[m] ],
  368.     DTFTtoZ[x_] :>
  369.         Transform[ x, 0, Infinity ],
  370.     DTFTtoZ[x_, w_, z_] :>
  371.         Transform[ substituteforw[x, w, -I Log[z]], 0, Infinity ],
  372.     myDTFT[ZTransData[x_, zrest__], rest__] :>
  373.         x,
  374.     substituteforw[x_, w_Symbol, w0_, opstr_, op_] :>
  375.         op[substituteforw[x, w, w - w0], substituteforw[x, w, w + w0]],
  376.     substituteforw[x_, w_Symbol, -w] :>
  377.         If [ Count[{x}, f_[p___][a___], Infinity] > 0,
  378.              Rev[w][posttransform[x]],          (* operators present *)
  379.              posttransform[x] /. w -> -w ] ,    (* no operators *)
  380.     substituteforw[x_, w_Symbol, wnew_] :>
  381.         ( posttransform[x] /. w -> wnew ),
  382.     upsamplemD[x_, {}, w_List, l_?MatrixQ] :>
  383.         ( posttransform[x] /. ReplaceWith[w, l . w] ),
  384.     upsamplemD[x_, nlist_, w_, l_] :>
  385.         upsamplemD[ posttransform[x], nlist, w, l ],
  386.     downsamplemD[x_, {}, w_, m_?MatrixQ] :>
  387.         postdownsamplemD[x, w, m],
  388.     downsamplemD[x_, nlist_, w_, m_] :>
  389.         downsamplemD[ posttransform[x], nlist, w, m ]
  390. }
  391.  
  392. ZtoDTFT[ SignalProcessing`Digital`ZTransform`Private`forwardz[x_, zrest__], n_, w_, z_, rest__ ] :=
  393.     DTFTtoZ[ myDTFT[x, n, w, z, rest], w, z ]
  394. ZtoDTFT[ x_, n_, w_, z_, rest__ ] := x
  395.  
  396.  
  397. (*  Printed forms of supported routines for Dialogue -> All option  *)
  398.  
  399. Format[ myDTFT[ x_, n_, w_, z_, options_, flag_ ] ] :=
  400.     SequenceForm[ ColumnForm[{"DTFT",
  401.                   "     " ~StringJoin~ ToString[n]}],
  402.               { x } ]
  403.  
  404. Format[ derivative[x_, w_Symbol, k_] ] := 
  405.     SequenceForm[ ColumnForm[{"D" Superscript[k],
  406.                   "  " ~StringJoin~ ToString[w]}],
  407.               { x } ]
  408.  
  409. Format[    substituteforw[x_, w_, wnew_] ] :=
  410.     SequenceForm[ {x}, Subscript[w -> wnew] ]
  411.  
  412. Format[    substituteforw[x_, w_, w0_, opstr_, op_] ] :=
  413.     SequenceForm[ {x}, Subscript[w -> w - w0], opstr,
  414.               {x}, Subscript[w -> w + w0] ]
  415.  
  416. Format[    DTFTtoZ[x_, w_, z_] ] :=
  417.     SequenceForm[ {x}, Subscript[w -> -I Log[z]] ]
  418.  
  419. Format[ downsample[x_, w_, m_] ] :=
  420.     SequenceForm[ {x}, Subscript[w -> ((w + 2 Pi "k")/m)] ]
  421.  
  422. Format[ downsamplemD[x_, n_, w_, m_] ] :=
  423.     SequenceForm[ {x}, Subscript["aliased by " m] ]
  424.  
  425. Format[ upsamplemD[x_, n_, w_, m_] ] :=
  426.     SequenceForm[ {x}, Subscript[w -> m . w] ]
  427.  
  428.  
  429. Format[ z ] := "z"
  430.  
  431.  
  432. (*  S T A T E  *)
  433.  
  434. initDTFTstate[] := True
  435. nullDTFTstate[] := False
  436.  
  437.  
  438. (*  B E G I N     D T F T     R U L E     B A S E  *)
  439.  
  440.  
  441. DTFTRules = {}
  442.  
  443.  
  444. OriginalDTFTRules = {
  445.  
  446.  
  447. (*  A.    Lookup rules        *)
  448. (*    1.  Pulse        *)
  449. myDTFT[ Pulse[L_, n_ + n0_.], n_, w_, z_, options_, flag_ ] :>
  450.     Periodic[2 Pi, w][Exp[I w n0] Exp[- I w (L - 1) / 2] Dirichlet[L, w]] /;
  451.     FreeQ[{L, n0}, n],
  452.  
  453. myDTFT[ A_. Step[n_ + n0_.] - A_. Step[n_ + n1_.] + x_, n_, w_, z_, options_, flag_ ] :>
  454.     myDTFT[ A Pulse[Abs[n1 - n0], n + Min[-n0, -n1]],
  455.         n, w, z, options, flag ] +
  456.     myDTFT[ x, n, w, z, options, flag ] /;
  457.     FreeQ[{n0, n1}, n],
  458.  
  459. (*    2.  Step        *)
  460. myDTFT[ Step[n_ + n0_.], n_, w_, z_, options_, flag_ ] :> 
  461.     Block [    {l = unique["l"]},
  462.         Exp[I w n0] / ( 1 - Exp[-I w] ) + 
  463.           Pi Exp[I w n0] Summation[l, -Infinity, Infinity, 1]
  464.                        [Delta[w + 2 Pi l]] ] /;
  465.     FreeQ[n0, n],
  466.  
  467. (*    3.  Summation forms    *)
  468. myDTFT[ Summation[k_, -Infinity, Infinity, 1][x_],
  469.         n_, w_, z_, options_, flag_ ] :>
  470.     Summation[k,-Infinity,Infinity,1][
  471.         myDTFT[x, n, w, z, options, flag] ] /;
  472.     FreeQ[k, n],
  473.  
  474. myDTFT[ Summation[k_, 0, L_, 1][a_. Exp[I m_ Pi k_ n_]],
  475.         n_, w_, z_, options_, flag_ ] :>
  476.     2 Pi Summation[k, -Infinity, Infinity,1][a Delta[w - 2 Pi k/(L+1)]] /;
  477.     FreeQ[{a, k, L}, n] && SameQ[m, 2/(L+1)],
  478.  
  479. (*    4.  Infinite extent cosine    *)
  480. myDTFT[ Cos[w0_. n_ + b_.], n_, w_, z_, options_, flag_ ] :> 
  481.     Block [    {l = unique["l"]},
  482.         Pi Summation[l,-Infinity,Infinity,1]
  483.                 [Exp[I b w / w0] Delta[w - w0 + 2 Pi l] +
  484.                  Exp[- I b w / w0] Delta[w + w0 + 2 Pi l]] ] /;
  485.     FreeQ[{b, w0}, n],
  486.  
  487. (*    5.  Infinite extent sine    *)
  488. myDTFT[ Sin[w0_. n_ + b_.], n_, w_, z_, options_, flag_ ] :> 
  489.     Block [    {l = unique["l"]},
  490.         Pi Summation[l,-Infinity,Infinity,1]
  491.                 [- I Exp[I b w / w0] Delta[w - w0 + 2 Pi l] +
  492.                  I Exp[I b w / w0] Delta[w + w0 + 2 Pi l]] ] /;
  493.     FreeQ[{b, w0}, n],
  494.  
  495. (*    6.  Constant function        *)
  496. myDTFT[ 1, n_, w_, z_, options_, flag_ ] :> 
  497.     Block [    {l = unique["l"]},
  498.         2 Pi Summation[l,-Infinity,Infinity,1][Delta[w + 2 Pi l]] ],
  499.  
  500. (*    7.  Sinc and other 1/n        *)
  501. myDTFT[ (-1)^(n_ + a_.) / ( n_ - 1 ), n_, w_, z_, options_, flag_ ] :>
  502.     Periodic[2 Pi, w][ (-1)^a w Exp[- I w] CPulse[2 Pi, w + Pi] ] /;
  503.     FreeQ[a, n],
  504.  
  505. myDTFT[ 1/n_, n_, w_, z_, options_, flag_ ] :>
  506.     Periodic[2 Pi, w][ - 2 Pi I CStep[w] ],
  507.  
  508. myDTFT[ Sinc[a_. n_ + b_.], n_, w_, z_, options_, flag_ ] :> 
  509.     Periodic[2 Pi, w][ Exp[I b/a] Pi CPulse[2 Abs[a], w + Abs[a]] /
  510.                Abs[a] ] /;
  511.     FreeQ[{a,b}, n],
  512.  
  513. myDTFT[ (Sinc[a_. n_ + b_.])^2, n_, w_, z_, options_, flag_ ] :> 
  514.     Periodic[2 Pi, w][ Exp[I b / a] (Pi / a)^2 ( Unit[-2][w + 2 a] -
  515.                2 Unit[-2][w] + Unit[-2][w - 2 a] ) ] /;
  516.     FreeQ[{a,b}, n],
  517.  
  518.  
  519. (*  B.    Properties        *)
  520. (*    1.  Multiplication by constant  *)
  521. myDTFT[ a_ x_, n_, w_, z_, options_, flag_ ] :>
  522.     a myDTFT[x, n, w, z, options, flag] /;
  523.     FreeQ[a, n],
  524.  
  525. (*    2.  Additivity property    *)
  526. myDTFT[ x_ + y_, n_, w_, z_, options_, flag_ ] :> 
  527.     myDTFT[x, n, w, z, options, flag] + myDTFT[y, n, w, z, options, flag],
  528.  
  529. (*    3.  Multiplication by n    *)
  530. myDTFT[ n_^a_. b_, n_, w_, z_, options_, flag_ ] :>
  531.     Block [ {dialogue},
  532.                 dialogue = SameQ[ Replace[Dialogue, options], All ];
  533.         Assuming[ a > 0, dialogue ];
  534.         If [ dialogue, Print[ "where ", a, " is an integer" ] ];
  535.         derivative[myDTFT[b, n, w, z, options, flag], w, a] ] /;
  536.     FreeQ[a,n] && Implies[ NumberQ[a], IntegerQ[a] && a > 0 ],
  537.  
  538. (*    4.  Delay by m        *)
  539. myDTFT[ f_. Step[n_ + m_], n_, w_, z_, options_, flag_ ] :>
  540.     Exp[I m w] myDTFT[(f /. n->(n - m)) Step[n], n, w, z, options, flag] /;
  541.     FreeQ[m, n],
  542.  
  543. (*    5.  Expand shifted sinusoids    *)
  544. myDTFT[ Cos[a_. n_ + b_] x_., n_, w_, z_, options_, flag_ ] :> 
  545.     myDTFT[ Cos[b] Cos[a n] x - Sin[b] Sin[a n] x,
  546.         n, w, z, options, flag ] /;
  547.     FreeQ[{a, b}, n],
  548.  
  549. myDTFT[ Sin[a_. n_ + b_] x_., n_, w_, z_, options_, flag_ ] :> 
  550.     myDTFT[ Sin[b] Cos[a n] x + Cos[b] Sin[a n] x,
  551.         n, w, z, options, flag ] /;
  552.     FreeQ[{a, b}, n],
  553.  
  554. (*    6.  Multiplication by complex exponential    *)
  555. myDTFT[ a_. Exp[Complex[0, w0_] w1_. n_ + b_.], n_, w_, z_, options_, flag_ ] :>
  556.     Exp[b] substituteforw[ myDTFT[a, n, w, z, options, flag],
  557.                    w, w - w0 w1 ] /;
  558.     FreeQ[{b, w0, w1}, n],
  559.  
  560. myDTFT[ a_. Cos[w0_. n_], n_, w_, z_, options_, flag_ ] :>
  561.     substituteforw[ myDTFT[a, n, w, z, options, flag] / 2,
  562.             w, w0, "+", Plus ] /;
  563.     FreeQ[w0, n],
  564.  
  565. myDTFT[ a_. Sin[w0_. n_], n_, w_, z_, options_, flag_ ] :>
  566.     substituteforw[ -I/2 myDTFT[a, n, w, z, options, flag],
  567.             w, w1, "-", Subtract ] /;
  568.     FreeQ[w0, n],
  569.  
  570. (*  C.    Structures        *)
  571. (*      -.  An operator independent of n---  take transform of f    *)
  572. myDTFT[ operator_[params__][f_], n_, w_, z_, options_, flag_ ] :>
  573.     operator[params][ myDTFT[f, n, w, z, options, flag] ] /;
  574.     FreeQ[{params}, n],
  575.  
  576. (*    1.  Upsampling        *)
  577. myDTFT[ Upsample[l_, n_] [ x_ ], n_, w_, z_, options_, flag_ ] :>
  578.     substituteforw[ myDTFT[x, n, w, z, options, flag], w, l w],
  579.  
  580. myDTFT[ Upsample[l_, nvars_List] [ x_ ], n_, w_, z_, options_, flag_ ] :>
  581.     Block [    {dim, expr, i},
  582.         dim = Length[nvars];
  583.         expr = x;
  584.         For [ i = 1, i <= dim, i++,
  585.               expr = Upsample[ l[[i,i]], nvars[[i]] ][expr] ];
  586.         myDTFT[ expr, n, w, z, options, flag ] ] /;
  587.     MemberQ[nvars, n] && diagonalResamplingMatrixQ[l],
  588.  
  589. myDTFT[ Upsample[l_, nvars_List] [ x_ ], n_, w_, z_, options_, flag_ ] :>
  590.     upsampleSetUp[l, nvars, x, n, w, z, options, flag ] /;
  591.     mDresampleCheck[Upsample, n, nvars, nList, l],
  592.  
  593. myDTFT[ upsamplemD[ x_, nlist_, wlist_, l_], n_, w_, z_, options_, flag_ ] :>
  594.     upsamplemD[ myDTFT[x, n, w, z, options, flag ],
  595.             If [ MemberQ[nlist, n], Complement[nlist, {n}], nlist ],
  596.             wlist, l ],
  597.  
  598. (*    2.  Downsampling    *)
  599. myDTFT[ Downsample[m_, n_] [ x_ ], n_, w_, z_, options_, flag_ ] :>
  600.     downsample[ myDTFT[x, n, w, z, options, flag], w, m ],
  601.  
  602. myDTFT[ Downsample[m_, nvars_List] [ x_ ], n_, w_, z_, options_, flag_ ] :>
  603.     Block [    {dim, expr, i},
  604.         dim = Length[nvars];
  605.         expr = x;
  606.         For [ i = 1, i <= dim, i++,
  607.               expr = Downsample[ m[[i,i]], nvars[[i]] ][expr] ];
  608.         myDTFT[expr, n, w, z, options, flag] ] /;
  609.     MemberQ[nvars, n] && diagonalResamplingMatrixQ[m],
  610.  
  611. myDTFT[ Downsample[m_, nvars_List] [ x_ ], n_, w_, z_, options_, flag_ ] :>
  612.     downsampleSetUp[m, nvars, x, n, w, z, options, flag ] /;
  613.     mDresampleCheck[Downsample, n, nvars, nList, m],
  614.  
  615. myDTFT[ downsamplemD[ x_, nlist_, wlist_, m_], n_, w_, z_, options_, flag_ ] :>
  616.     downsamplemD[ myDTFT[x, n, w, z, options, flag ],
  617.               If [ MemberQ[nlist, n], Complement[nlist, {n}], nlist ],
  618.               wlist, m ],
  619.  
  620. (*  D.    Strategies        *)
  621. (*    1.  Multiply out terms in cosine, sin, and exponential functions.  *)
  622. myDTFT[ Exp[a_] x_., n_, w_, z_, options_, flag_ ] :>
  623.     myDTFT[ Exp[Expand[a]] x, n, w, options, flag ] /;
  624.     ! SameQ[Exp[a], Exp[Expand[a]]],
  625.  
  626. myDTFT[ Cos[a_] x_., n_, w_, z_, options_, flag_ ] :>
  627.     myDTFT[ Cos[Expand[a]] x, n, w, options, flag ] /;
  628.     ! SameQ[Cos[a], Cos[Expand[a]]],
  629.  
  630. myDTFT[ Sin[a_] x_., n_, w_, z_, options_, flag_ ] :>
  631.     myDTFT[ Sin[Expand[a]] x, n, w, options, flag ] /;
  632.     ! SameQ[Sin[a], Sin[Expand[a]]],
  633.  
  634. (*    2.  Expand out products like (n + 1)(n + 2) ...    *)
  635. myDTFT[ x_, n_, w_, z_, options_, flag_ ] :>
  636.     myDTFT[ Distribute[x], n, w, z, options, flag ] /;
  637.     SameQ[Head[x], Times] && ! SameQ[x, Distribute[x]],
  638.  
  639. (*    3.  Use the z-transform    *)
  640. myDTFT[ x_, n_, w_, z_, options_, True ] :> 
  641.     Block [    {myz, nullstate, ztrans},
  642.         nullstate = nullDTFTstate[];
  643.         myz = SignalProcessing`Digital`ZTransform`Private`MyZByPass;
  644.         ztrans = myz[ x, n, z, nList, zList,
  645.                   {Dialogue -> False} ~Join~ options ];
  646.         If [ InformUserQ[options],
  647.              Print[ "( after taking the z-transform of" ];
  648.              Print[ "  ", x ];
  649.              Print[ "  which yields" ];
  650.              Print[ "  ", ztrans ];
  651.              Print[ "  and adjusting the result )" ] ];
  652.         If [ SameQ[Head[ztrans], ZTransData],
  653.              ( TheFunction[ztrans] /. z -> Exp[I w] ),
  654.              myDTFT[ MapAll[ ZtoDTFT[#1, n, w, z, options, nullstate]&,
  655.                      ztrans ],
  656.                  n, w, z, options, nullstate ] ] ]
  657.  
  658. }
  659.  
  660. (*  E N D     R U L E     B A S E  *)
  661.  
  662. (*  E N D     F O R W A R D     D T F T  *)
  663.  
  664.  
  665. (*  B E G I N     I N V E R S E     D T F T  *)
  666.  
  667. InvDTFTransform::noncont =
  668.     "Discrete operators and functions were converted into continuous ones."
  669.  
  670. (*  Conversion from operator notation to function notation  *)
  671. Unprotect[InvDTFT]
  672. InvDTFT/: TheFunction [ InvDTFT[w_, n_][f_] ] := InvDTFTransform[f, w, n]
  673. Protect[InvDTFT]
  674.  
  675. (*  Give it options  *)
  676. InvDTFTransform/: Options[InvDTFTransform] :=
  677.     displayFixUp[ invFixUp,
  678.               Join[ { Definition -> False,
  679.                   Dialogue -> False,
  680.                   Terms -> False },
  681.                 Options[InvZTransform] ] ]
  682.  
  683. (*  Input a function and test to see if needs to be sent to the    *)
  684. (*  multidimensional or one-dimensional inverse Fourier transform. *)
  685. InvDTFTransform[ DTFTData[ trans_, FVariables[vars_] ] ] :=
  686.     InvDTFTransform [ trans, vars ]
  687.  
  688. InvDTFTransform[ f_ ] :=
  689.     Message[ Transform::novariables, "w (frequency)", GetVariables[f] ]
  690.  
  691. InvDTFTransform[ f_, w_ ] :=
  692.     InvDTFTransform[ f, w, DummyVariables[Length[w], Global`n] ]
  693.  
  694. InvDTFTransform[ f_, w_, n_, options___ ] :=
  695.     Message[ Transform::badvar, "continuous-frequency",
  696.          InvDTFTransform, w ] /;
  697.     ! validvarQ[w]
  698.  
  699. InvDTFTransform[ f_, w_, n_, options___ ] :=
  700.     Message[ Transform::badvar, "discrete-time", InvDTFTransform, n ] /;
  701.     ! validvarQ[n]
  702.  
  703. InvDTFTransform[ f_, varin_, varout_, options___ ] :=
  704.     Block [    {freqfun, invfun, op, vars, z, zvars},
  705.         vars = If [ Length[varin] != Length[varout],
  706.                 DummyVariables[Length[varin], Global`n],
  707.                 varout ];
  708.  
  709.         zvars = DummyVariables[Length[varin], z];
  710.  
  711.         freqfun = ToContinuous[f];
  712.         If [ ! SameQ[freqfun, f],
  713.              Message[InvDTFTransform::noncont] ];
  714.  
  715.         wList = ToList[varin];
  716.         nList = ToList[vars];
  717.         zList = ToList[zvars];
  718.  
  719.         op = ToList[options] ~Join~ Options[InvDTFTransform];
  720.         multidFlag = ( Length[varin] > 1 );
  721.         invfun = If [ multidFlag, 
  722.                   multiDInvDTFT[freqfun, varin, vars, zvars, op],
  723.                   MyInvDTFT[freqfun, varin, vars, zvars, op] ];
  724.  
  725.         If [ InformUserQ[op],
  726.              Scan[invexplain, invfun, Infinity] ];
  727.  
  728.         invfun ]
  729.  
  730. (*  Multidimensional inverse discrete-time Fourier transform  *)
  731. multiDInvDTFT[ f_, wl_, nl_, zl_, oplist_ ] :=
  732.     Block [    {fun = f, i, len, newop},
  733.         len = Length[wl];
  734.         newop = MultiDLookup[oplist, wl, nl];
  735.         For [ i = 1, i <= len, i++,
  736.               fun = MyInvDTFT[fun, wl[[i]], nl[[i]], zl[[i]], newop] ];
  737.         fun ]                
  738.  
  739. (*  One-dimensional inverse discrete-time Fourier transform  *)
  740. MyInvDTFT[ f_, w_, n_, z_, op_ ] :=
  741.     Block [ {newexpr = myinvDTFT[f, w, n, z, invFixUp[op], {True, True}],
  742.          newinvdtftrules, oldexpr = Null, trace},
  743.  
  744.         (* generate the inverse DTFT rules *)
  745.         newinvdtftrules = TransformFixUp[ w, n, op, myinvDTFT, False,
  746.                           InvDTFTransform, Null, Null ];
  747.         InvDTFTRules = Join[ newinvdtftrules, OriginalInvDTFTRules ];
  748.  
  749.         (* determine the level of dialogue *)
  750.         trace = SameQ[ Replace[Dialogue, op], All ];
  751.  
  752.         (* take the 1-D inverse DTFT *)
  753.         While [ ! SameQ[newexpr, oldexpr], 
  754.             If [ trace, Print[newexpr]; Print[ "which becomes" ] ];
  755.             oldexpr = newexpr;
  756.             newexpr = oldexpr /. ( myinvDTFT[a__] :>
  757.                     Replace[ myinvDTFT[a], InvDTFTRules ] ) ];
  758.         If [ trace, Print[newexpr] ];
  759.  
  760.         (* fix up the result *)
  761.         newexpr = postinvtransform[oldexpr];
  762.         While [ ! SameQ[newexpr, oldexpr], 
  763.             If [ trace, Print[newexpr]; Print[ "which becomes" ] ];
  764.             oldexpr = newexpr;
  765.             newexpr = postinvtransform[oldexpr] ];
  766.  
  767.         (* simplify expression if requested *)
  768.         If [ Replace[Simplify, op],
  769.              newexpr = SPSimplify[newexpr,
  770.                       Dialogue -> trace,
  771.                       Variables -> nList] ];
  772.  
  773.         newexpr ]
  774.  
  775.  
  776. (*  Supporting Routines  *)
  777.  
  778. breakItUp[ a_. f_, w_, k_ ] := { a, f } /;
  779.     FreeQ[a, k] && ! FreeQ[f, k]
  780.  
  781. checkFactoring[ a_, w_, k_ ] := breakItUp[ Factor[a], w, k ]
  782.  
  783. checkFactoringQ[ a_, w_, k_ ] :=
  784.     ! SameQ[ First[ checkFactoring[ a, w, k ] ], 1 ]
  785.  
  786. candidateList = {}
  787. periodicX = {}
  788. nonPeriodicX = {}
  789.  
  790. downShiftPair[x_, y_, m_] := 
  791.     Block [    {diff, ratio},
  792.         diff = x - y;
  793.         If [ Apply[SameQ, diff],
  794.              ratio = First[diff] / (2 Pi);
  795.              IntegerQ[ratio] && ( -m <= ratio <= m ) ] ] /;
  796.     Length[x] == Length[y]
  797.  
  798. downsampled1DcheckQ[ x_Plus, w_, m_ ] :=
  799.     Block [    {i, j, length, lengthList, newx, rule},
  800.         newx = MyCollectAll[x, w];
  801.         periodicX = {};
  802.         rule = a_. Periodic[2 Pi m, w][f_] :>
  803.                Block [ {}, AppendTo[periodicX, a f]; 0 ];
  804.         nonPeriodicX = Map[ Replace[#, rule]&, newx ];
  805.  
  806.         (*  Determine which functions are shifted versions of  *)
  807.         (*  one another--- first compute all shift factors for *)
  808.         (*  each candidate and then compare them pairwise      *)
  809.         shiftFactors =
  810.             Expand[ Map[ GetAllShiftFactors[#, w]&, periodicX ] ];
  811.         length = Length[ shiftFactors ];
  812.         lengthList = Map[ Length, shiftFactors ];
  813.         For [ i = 1, i <= length, i++,
  814.               candidateList = {i};
  815.               For [ j = 1, j <= length, j++,
  816.                 If [ (i != j) &&
  817.                    downShiftPair[ shiftFactors[[i]],
  818.                           shiftFactors[[j]], m ],
  819.                  AppendTo[candidateList, j] ] ];
  820.               If [ SameQ[Length[candidateList], m], Break[] ] ];
  821.         SameQ[ Length[candidateList], m ] ]
  822.  
  823. invexplain[ myinvDTFT[ x_, w_, rest__ ] ] :=
  824.     Message[ Transform::incomplete,
  825.          "inverse discrete-time Fourier transform", x, w ]
  826.  
  827. (*  invFixUp --  removes redundancy and TransformLookup in options list  *)
  828. invFixUp[ op_ ] := { Definition -> Replace[Definition, op],
  829.              Dialogue -> Replace[Dialogue, op],
  830.              Terms -> Replace[Terms, op] }
  831.  
  832. (*  invertDownsampledDTFT  *)
  833. invertDownsampledDTFT[ x_, w_, n_, z_, options_, flag_, m_ ] :=
  834.     Block [ {},
  835.         If [ ! IntegerQ[m],
  836.                      If [ SameQ[ Replace[Dialogue, options], All ],
  837.                   Print[ "assuming ", m, " is an integer" ] ] ];
  838.         Abs[m] Downsample[m, n][
  839.              myinvDTFT[x /. w -> m w, w, n, z, options, flag] ] ] /;
  840.     Implies[ NumberQ[m], IntegerQ[m] ]
  841.  
  842. (*  invertResampledDTFT  *)
  843. invertResampledDTFT[ x_, w_, n_, z_, options_, flag_, y_, f_ ] :=
  844.     invertResampledDTFT[ x, w, n, z, options, flag,
  845.                  y, Numerator[f], Denominator[f] ]
  846.  
  847. invertResampledDTFT[ x_, w_, n_, z_, options_, flag_, y_, m_, l_ ] :=
  848.     Block [    {extraterms = 0, newx, upfact},
  849.         upfact = If [ SameQ[m, 2], l, 2 l ];
  850.         If [ IntegerQ[upfact] && upfact < 0, upfact = -upfact ];
  851.         If [ ! IntegerQ[upfact],
  852.                      If [ SameQ[ Replace[Dialogue, options], All ],
  853.                   Print[ "assuming ", upfact, " is an integer" ] ] ];
  854.         newx = x /. w -> (w/upfact);
  855.         If [ SameQ[y, 0],
  856.              extraterms = myinvDTFT[y, w, n, z, options, flag] ];
  857.         extraterms +
  858.         Upsample[upfact, n][
  859.             myinvDTFT[newx, w, n, z, options, flag] ] ] /;
  860.     ( SameQ[m, 2] || SameQ[m, 1] ) && FreeQ[l, w] && (! SameQ[l, 1]) &&
  861.       Implies[ NumberQ[l], IntegerQ[l] ]
  862.  
  863. invertResampledDTFT[ x_, w_, n_, z_, options_, flag_, y_, m_Integer, 1 ] :=
  864.     Block [    {extraterms},
  865.         extraterms = nonPeriodicX;
  866.         If [ Length[candidateList] < m/2,
  867.              indices = Intersection[ Range[m/2], candidateList ];
  868.              extraterms += Apply[ Plus, periodicX[[indices]] ] ];
  869.         thefun = periodicX[[ First[candidateList] ]];
  870.         If [ ! SameQ[extraterms, 0],
  871.              extraterms =
  872.             myinvDTFT[extraterms, w, n, z, options, flag] ];
  873.         invertDownsampledDTFT[ thefun, w, n, z, options, flag, m/2 ] +
  874.           extraterms ] /;
  875.     EvenQ[m] && ( m > 2 ) &&
  876.       downsampled1DcheckQ[ Periodic[m Pi, w][x] + y, w, m/2 ]
  877.  
  878. invertResampledDTFT[ x_, w_, n_, z_, options_, flag_, y_, m_Integer, 1 ] :=
  879.     invertResampledDTFT[x, w, n, z, options, flag, y, 2 m, 2 ] /;
  880.     OddQ[m]
  881.  
  882.  
  883. postinvtransform[x_] := ReplaceRepeated[x, postinvrules]
  884.  
  885. postinvrules = {
  886.  
  887.     (*    4.  Apply the definition        *)
  888.     myinvDTFT[ a_ Pulse[l_, w_ + s_.], w_, n_, z_, options_, flag_ ] :>
  889.         TransformDialogue[ Definition, options, Integrate,
  890.             a Pulse[l, w + s],
  891.             Integrate[ a Exp[I w n], {w, -s, l - s}] / (2 Pi) ] /;
  892.         Replace[Definition, options] && N[(Pi > s) && (Pi < l - s)],
  893.  
  894.     myinvDTFT[ x_, w_, n_, z_, options_, flag_ ] :>
  895.         TransformDialogue[ Definition, options, Integrate, x,
  896.                Integrate[x Exp[I w n], {w, -Pi, Pi}] / (2 Pi) ] /;
  897.         Replace[Definition, options],
  898.  
  899.     substituteforn[x_, n_Symbol, -n_] :>
  900.         If [ Count[{x}, f_[p___][a___], Infinity] > 0,
  901.              Rev[n][posttransform[x]],          (* operators present *)
  902.              posttransform[x] /. n -> -n ] ,    (* no operators *)
  903.  
  904.     substituteforn[x_, n_, newn_] :> ( x /. n -> newn )
  905. }
  906.  
  907. (* Routines to detect upsampling in the frequency response  *)
  908.  
  909. upsamplefactorQ[ f_ ] := (! SameQ[f, 0]) && (! SameQ[f, 1]) && (! SameQ[f, -1])
  910. upsampledQ[ x_, w_ ] := upsamplefactorQ[ UpsampleFactor[x, w, GetAllFactors] ]
  911.  
  912. upsampledmDQ[ x_, w_ ] :=
  913.     Block [    {cond},
  914.         cond = ! MyFreeQ[ GetAllShiftFactors[x, w], wList ];
  915.         If [ cond,
  916.              Needs[ "SignalProcessing`Digital`MDDTFT`" ];
  917.              cond = upsampleSystemQ[ x, wList ] ];
  918.         cond ]
  919.  
  920.  
  921. (*  Printed forms of supported routines for Dialogue -> All option  *)
  922.  
  923. Format[ myinvDTFT[ x_, w_, n_, z_, options_, flag_ ] ] :=
  924.     SequenceForm[ ColumnForm[{"DTFT" Superscript[-1],
  925.                   "     " ~StringJoin~ ToString[w]}],
  926.               { x } ]
  927.  
  928. Format[    substituteforn[x_, n_, newn_] ] :=
  929.     SequenceForm[ {x}, Subscript[n -> newn] ]
  930.  
  931.  
  932. (*  B E G I N     R U L E     B A S E  *)
  933.  
  934.  
  935. InvDTFTRules = {}
  936.  
  937.  
  938. OriginalInvDTFTRules = {
  939.  
  940.  
  941. myinvDTFT[ Periodic[2 Pi, w_][x_], w_, n_, z_, options_, flag_ ] :>
  942.     myinvDTFT[x, w, n, z, options, flag],
  943.  
  944.  
  945. (*  A.    Lookup Rules        *)
  946. (*    1.  Summation forms    *)
  947. myinvDTFT[ Summation[k_,-Infinity,Infinity,1][a_. Delta[w_+w0_.+b_. Pi k_] +
  948.                     a_. Delta[w_ + A_. + c_. Pi k_] + x_.],
  949.        w_, n_, z_, options_, flag_ ] :>
  950.     a Cos[A n] / Pi + Summation[k, -Infinity, Infinity, 1][x] /;
  951.     ( A == - w0 ) && ( SameQ[b, 2] || SameQ[b, -2] ) &&
  952.       ( SameQ[c, 2] || SameQ[c, -2] ) && FreeQ[{a, k, w0}, w],
  953.  
  954. myinvDTFT[ Summation[k_,-Infinity,Infinity,1][a_. Delta[w_+w0_.+b_. Pi k_] +
  955.                       d_. Delta[w_+A_.+c_. Pi k_] + x_.],
  956.        w_, n_, z_, options_, flag_ ] :>
  957.     (a I) Sin[A n] / Pi + Summation[k, -Infinity, Infinity, 1][x] /;
  958.     ( A == -w0 ) && ( a == -d ) && ( SameQ[b, 2] || SameQ[b, -2] ) &&
  959.       ( SameQ[c, 2] || SameQ[c, -2] ) && FreeQ[{a, A, w0}, w],
  960.  
  961. myinvDTFT[ Summation[k_, -Infinity, Infinity, 1][
  962.                     a_. Delta[w_+w0_.+b_. Pi k_] + x_.],
  963.        w_, n_, z_, options_, flag_ ] :>
  964.     a Exp[- I w0 n] / (2 Pi) + Summation[k, -Infinity, Infinity, 1][x] /;
  965.     ( SameQ[b, 2] || SameQ[b, -2] ) && FreeQ[{a, k, w0}, w],
  966.  
  967. myinvDTFT[ Summation[k_,-Infinity,Infinity,1][a_. Delta[w_ + b_. Pi k_ / L_]],
  968.         w_, n_, z_, options_, flag_ ] :>
  969.     L/(2 Pi) Summation[k,-Infinity,Infinity,1][a Exp[I 2 Pi k n / (L+1)]] /;
  970.     FreeQ[{a, k, L}, w] && ( SameQ[b, -2] || SameQ[b, 2] ),
  971.  
  972. (*    2.  Constant function    *)
  973. myinvDTFT[ a_, w_, n_, z_, options_, flag_ ] :> a Impulse[n] /; FreeQ[a, w],
  974.  
  975. (*    3.  Pulse in continuous frequency    *)
  976. myinvDTFT[ CPulse[L_, a_. w_ + w1_.], w_, n_, z_, options_, flag_ ] :>
  977.     L/a Exp[I w1 n / a] Exp[- I L n / (2 a)] Sinc[L n / (2 a)] / (2 Pi) /;
  978.     FreeQ[{a, L, w1}, w],
  979.  
  980. myinvDTFT[ A_. CStep[w_ + w0_.] - A_. CStep[w_ + w1_.] + x_.,
  981.        w_, n_, z_, options_, flag_ ] :>
  982.     myinvDTFT[ A CPulse[w0 - w1, w + w0], w, n, z, options, flag ] +
  983.       myinvDTFT[ x, w, n, z, options, flag ] /;
  984.     FreeQ[{w0, w1}, w],
  985.  
  986. (*    4.  Step function    *)
  987. myinvDTFT[ 1 / ( 1 - Exp[- I w_] ) +
  988.        p1_. Summation[k_,-Infinity,Infinity,1][p2_. Delta[w_ + b_. Pi k_]],
  989.        w_, n_, z_, options_, flag_] :> 
  990.     Step[n] /;
  991.     FreeQ[k, w] && SameQ[p1 p2, Pi] && ( SameQ[b, 2] || SameQ[b, -2] ),
  992.  
  993. (*     5.  Dirichlet kernel (aliased sinc)    *)
  994. myinvDTFT[ Dirichlet[L_, w_ + w0_.], w_, n_, z_, options_, flag_ ] :>
  995.     Exp[- I w0 n] Pulse[L, n + (L - 1)/2] /;
  996.     FreeQ[{L, w0}, w],
  997.  
  998. (*    6.  Differentiator (high pass) filter    *)
  999. myinvDTFT[ w_ CPulse[2 Pi, w_ + Pi], w_, n_, z_, options_, flag_ ] :>
  1000.     -2 I Pi (-1)^n / n,
  1001.  
  1002.  
  1003. (*  B.    Properties                *)
  1004. (*      1.  Resampling---  all Periodic[2 Pi, w] operators have already  *)
  1005. (*               been processed by the rule base.         *)
  1006. myinvDTFT[ Summation[k_, 0, Lminus1_, 1][ Periodic[m_. Pi, w_][ x_ ] ],
  1007.        w_, n_, z_, options_, flag_ ] :>
  1008.     invertDownsampledDTFT[ x /. k -> 0, w, n, z, options, flag, m/2 ] /;
  1009.     ( SameQ[m/2, Lminus1 + 1] || SameQ[Abs[m/2], Lminus1 + 1] ) &&
  1010.     ( Count[ x, w/(m/2) + 2 Pi k/(m/2) + h_., Infinity ] > 0 ||
  1011.       Count[ x, (w + 2 Pi k + h_.)/(m/2), Infinity ] > 0 ),
  1012.  
  1013. myinvDTFT[ Periodic[m_. Pi, w_][ x_ ] + y_., w_, n_, z_, options_, flag_ ] :>
  1014.     invertResampledDTFT[ x, w, n, z, options, flag, y, Rationalize[m] ],
  1015.  
  1016. (*    2.  Multiplication by a constant    *)
  1017. myinvDTFT[ a_ x_, w_, n_, z_, options_, flag_ ] :>
  1018.     a myinvDTFT[x, w, n, z, options, flag] /;
  1019.     FreeQ[a, w],
  1020.  
  1021. (*    3.  Additivity                *)
  1022. myinvDTFT[ a_ + b_, w_, n_, z_, options_, flag_ ] :>
  1023.     myinvDTFT[a, w, n, z, options, flag] +
  1024.     myinvDTFT[b, w, n, z, options, flag],
  1025.  
  1026. (*    4.  Multiplication by a complex exp.    *)
  1027. myinvDTFT[ Exp[Complex[0, n1_] n0_. w_ + b_.] x_.,
  1028.        w_, n_, z_, options_, flag_ ] :>
  1029.     Exp[b] substituteforn[ myinvDTFT[x, w, n, z, options, flag],
  1030.                    n, n + n0 n1 ] /;
  1031.     FreeQ[{b, n0, n1}, w],
  1032.  
  1033. (*    4.  Summation                *)
  1034. myinvDTFT[ x_. / ( 1 - Exp[- I w_] ) +
  1035.        Pi Summation[k_,-Infinity,Infinity,1][Delta[w_ + b_. Pi k_]],
  1036.        w_, n_, z_, options_, flag_ ] :>
  1037.     Summation[k, -Infinity, Infinity, 1][
  1038.       myinvDTFT[x, w, n, z, options, flag] ] /;
  1039.     FreeQ[k, w] && ( SameQ[b, -2] || SameQ[b, 2] ),
  1040.  
  1041. myinvDTFT[ b_. Summation[k_,-Infinity,Infinity,1][a_],
  1042.        w_, n_, z_, options_, flag_ ] :>
  1043.     Block [    {f, newdtft, rest},
  1044.         {f, rest} = checkFactoring[a, w, k];
  1045.         newdtft = b f Summation[k,-Infinity,Infinity,1][Expand[rest]];
  1046.         myinvDTFT[ newdtft, w, n, z, options, flag ] ] /;
  1047.     checkFactoringQ[a, w, k],
  1048.  
  1049. (*    5.  Multdimensional upsampling    *)
  1050. myinvDTFT[ x_, w_, n_, z_, options_, flag_ ] :>
  1051.     undoUpsampling[ x, w, n, z, options, flag ] /;
  1052.     multidFlag && upsampledmDQ[ x, w ],
  1053.  
  1054. (*      6.  An operator independent of w---  take transform of f    *)
  1055. myinvDTFT[ operator_[params__][f_], w_, n_, z_, options_, flag_ ] :>
  1056.     operator[params][ myinvDTFT[f, w, n, z, options, flag] ] /;
  1057.     FreeQ[{params}, w],
  1058.  
  1059.  
  1060. (*  C.    Strategies                    *)
  1061. (*    1.  Collect all subexpressions in terms of w    *)
  1062. myinvDTFT[ x_, w_, n_, z_, options_, {True, flag2_} ] :>
  1063.     myinvDTFT[ MyCollectAll[x, w], w, n, z, options, { False, flag2 } ],
  1064.  
  1065. (*    2.  Expand out products like (w + 1)(w + 2) ...    *)
  1066. myinvDTFT[ x_, w_, n_, z_, options_, flag_ ] :>
  1067.     myDTFT[ Distribute[x], n, w, z, options, flag ] /;
  1068.     SameQ[Head[x], Times] && ! SameQ[x, Distribute[x]],
  1069.  
  1070. (*    3.  Take inverse z-transform            *)
  1071. myinvDTFT[ x_, w_, n_, z_, options_, {flag1_, True} ] :>
  1072.     Block [    {myinvz, newx, invztrans},
  1073.         newx = SPSimplify[ ( x /. w -> ( - I Log[z] ) ) ];
  1074.         myinvz =
  1075.           SignalProcessing`Digital`InvZTransform`Private`MyInvZByPass;
  1076.         invztrans = myinvz[ newx, z, n, zList, nList,
  1077.                     {Dialogue -> False} ~Join~ options ];
  1078.         If [ InformUserQ[options],
  1079.              Print[ "( after taking the inverse z-transform of" ];
  1080.              Print[ "  ", x ];
  1081.              Print[ "  which yields" ];
  1082.              Print[ "  ", invztrans ];
  1083.              Print[ "  and adjusting the result )" ] ];
  1084.         If [ InvalidInvZTransformQ[invztrans],
  1085.              myinvDTFT[ invztrans, w, n, z, options, {flag1, False} ],
  1086.              invztrans ] ],
  1087.  
  1088. (*        Clean up expression returned from inverse z-transform    *)
  1089. myinvDTFT[SignalProcessing`Digital`InvZTransform`Private`myinvz[x_, z_, rest__], w_, n_, z_, options_, flag_ ] :>
  1090.     Block [    {dtftform},
  1091.         dtftform = SPSimplify[ ( x /. z -> Exp[I w] ) ];
  1092.         myinvDTFT[dtftform, w, n, z, options, flag] ]
  1093.  
  1094. }
  1095.  
  1096.  
  1097. (*  E N D     R U L E     B A S E  *)
  1098.  
  1099. (*  E N D     I N V E R S E     D T F T  *)
  1100.  
  1101.  
  1102. (*  E N D     P A C K A G E  *)
  1103.  
  1104. End[]
  1105. EndPackage[]
  1106.  
  1107. If [ TrueQ[ $VersionNumber >= 2.0 ],
  1108.      On[ General::spell ];
  1109.      On[ General::spell1 ] ];
  1110.  
  1111.  
  1112. (*  H E L P     I N F O R M A T I O N  *)
  1113.  
  1114. Combine [ SPfunctions, { DTFTransform, InvDTFTransform } ]
  1115. Protect [ DTFTransform, InvDTFTransform ]
  1116.  
  1117.  
  1118. (*  E N D I N G     M E S S A G E  *)
  1119.  
  1120. Print["The discrete-time Fourier transform (DTFT) rule bases are loaded."]
  1121. Null
  1122.